survival hw 6

a

library(survival)
library(ClinicalTrialSummary)
data(ggas)
# group 0: Chemotherapy only; 1: Chemotherapy plus radiotherapy
dat <- ggas
dat$time <- dat$time * 365

# 暴露组： Chemotherapy only 对照组：Chemotherapy plus radiotherapy
dat$therapy <- ifelse(dat$group == 0, 1, 0) 

model_a <- coxph(Surv(time, event) ~ therapy, dat)
summary(model_a)

Call:
coxph(formula = Surv(time, event) ~ therapy, data = dat)

  n= 90, number of events= 82 

           coef exp(coef) se(coef)      z Pr(>|z|)
therapy -0.1051    0.9002   0.2233 -0.471    0.638

        exp(coef) exp(-coef) lower .95 upper .95
therapy    0.9002      1.111    0.5811     1.395

Concordance= 0.562  (se = 0.031 )
Likelihood ratio test= 0.22  on 1 df,   p=0.6
Wald test            = 0.22  on 1 df,   p=0.6
Score (logrank) test = 0.22  on 1 df,   p=0.6

结果解读

Cox模型结果显示，治疗（therapy）对事件风险的影响不显著（HR = 0.90, 95% CI: 0.58–1.40, p = 0.638），风险比接近1，说明治疗效果较弱且无统计学意义。模型的一致性指数为0.562，预测能力较低，三个检验（似然比、Wald和对数秩）p值均为0.6，进一步支持治疗的作用不显著。需更大样本或更明确的治疗方案验证其临床意义。

b

# 这里t由coxph直接传入了
model_b <- coxph(Surv(time, event) ~ therapy + tt(therapy), dat,
               tt = function(x, t, ...)x*log(t))
summary(model_b)

Call:
coxph(formula = Surv(time, event) ~ therapy + tt(therapy), data = dat, 
    tt = function(x, t, ...) x * log(t))

  n= 90, number of events= 82 

                coef exp(coef) se(coef)      z Pr(>|z|)   
therapy     -3.86619   0.02094  1.47707 -2.617  0.00886 **
tt(therapy)  0.64677   1.90935  0.24793  2.609  0.00909 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

            exp(coef) exp(-coef) lower .95 upper .95
therapy       0.02094    47.7602  0.001158    0.3786
tt(therapy)   1.90935     0.5237  1.174488    3.1040

Concordance= 0.576  (se = 0.03 )
Likelihood ratio test= 8.54  on 2 df,   p=0.01
Wald test            = 6.87  on 2 df,   p=0.03
Score (logrank) test = 7.59  on 2 df,   p=0.02

# KK书上引入时依变量的方法
dat.cp <- survSplit(dat,cut = dat$time[dat$event==1],
                    end = 'time',event = 'event',start = 'start',id='id')

dat.cp$logttherapy <- dat.cp$therapy *log(dat.cp$time)

coxph(Surv(dat.cp$start,dat.cp$time,dat.cp$event)~
        logttherapy+therapy+cluster(id),data=dat.cp)

Call:
coxph(formula = Surv(dat.cp$start, dat.cp$time, dat.cp$event) ~ 
    logttherapy + therapy, data = dat.cp, cluster = id)

                coef exp(coef) se(coef) robust se      z     p
logttherapy  0.64677   1.90935  0.24793   0.41179  1.571 0.116
therapy     -3.86619   0.02094  1.47707   2.47003 -1.565 0.118

Likelihood ratio test=8.54  on 2 df, p=0.01396
n= 3960, number of events= 82 

c

# KK书上关于时依变量的代码
# load("/Users/hcy/Zotero/storage/CE9IAU3S/addicts.rda")
# head(addicts)
#  
#  
# addicts.cp <- survSplit(addicts,cut = addicts$survt[addicts$status==1],end = 'survt',
#                        event = 'status',start = 'start')
#  
# head(addicts.cp)
#  
# addicts.cp$logtdose <- addicts.cp$dose*log(addicts.cp$survt)
#  
# addicts.cp[addicts.cp$id==1,]
# addicts[addicts$id==1,]
#  
# coxph(Surv(addicts.cp$start,addicts.cp$survt,addicts.cp$status)~
#         prison+dose+clinic+logtdose+cluster(id),data=addicts.cp)
# 
# addicts.cp365 <- survSplit(addicts,cut = 365,end = 'survt',
#                         event = 'status',start = 'start')
# 
# addicts.cp365$hv1 <- addicts.cp365$clinic*(addicts.cp365$start<365)
# 
# addicts.cp365$hv2 <- addicts.cp365$clinic*(addicts.cp365$start>=365)
# 
# Y365=Surv(addicts.cp365$start,addicts.cp365$survt, addicts.cp365$status)
# 
# coxph(Y365 ~ prison + dose + hv1 + hv2 + cluster(id), data=addicts.cp365)
# 
# coxph(Y365 ~ prison + dose + hv1 + hv2 + cluster(id), data=addicts.cp365,method = 'breslow')
# 
# coxph(Y365 ~ prison + dose + clinc + hv1 + hv2 + cluster(id), data=addicts.cp365,method = 'breslow')

tao <- unique(sort(dat$time[dat$event==1]))
l <- numeric(length(tao))
for (i in 1:length(tao)) {
  temp <- survSplit(Surv(time, event) ~ group, data= dat, cut=c(tao[i]),
                    episode= "tgroup", id="id",start='tstart')
  fit <- coxph(Surv(tstart, time, event) ~ group +
                 group:factor(tgroup), data=temp)
  l[i] <- logLik(fit)
}

Warning in agreg.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 1,2 ; beta may be infinite.

Warning in agreg.fit(X, Y, istrat, offset, init, control, weights = weights, :
Loglik converged before variable 2 ; beta may be infinite.

(taoHat <- tao[which.max(l)])

[1] 254

d

temp_d <- survSplit(Surv(time, event) ~ group, data= dat, cut=254,
                  episode= "tgroup", id="id")
temp_d$group1 <- ifelse(temp_d$group==0 & temp_d$tgroup==1,1,0)
temp_d$group2 <- ifelse(temp_d$group==0 & temp_d$tgroup==2,1,0)
fit_d2 <- coxph(Surv(tstart, time, event) ~ group1 + group2, data=temp_d)
summary(fit_d2)

Call:
coxph(formula = Surv(tstart, time, event) ~ group1 + group2, 
    data = temp_d)

  n= 150, number of events= 82 

          coef exp(coef) se(coef)      z Pr(>|z|)   
group1 -1.4229    0.2410   0.4327 -3.288  0.00101 **
group2  0.6423    1.9008   0.3048  2.107  0.03508 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

       exp(coef) exp(-coef) lower .95 upper .95
group1     0.241     4.1491    0.1032    0.5628
group2     1.901     0.5261    1.0460    3.4542

Concordance= 0.62  (se = 0.029 )
Likelihood ratio test= 17.79  on 2 df,   p=1e-04
Wald test            = 15.25  on 2 df,   p=5e-04
Score (logrank) test = 17.31  on 2 df,   p=2e-04